#include "fermiqcd.h"
// #include "fermiqcd_fermilab_action.h"
#include "fermiqcd_fermilab_action_sse.h"

void main(int argc, char** argv) {
  string info=argv[1];
  mpi.open_wormholes(argc, argv);
  define_base_matrices("FERMILAB");

  int box[4];
  coefficients coeff;

  string gauge_filename=          prompt(info, "INPUT_GAUGE_FILENAME", "COLD");
  string fermi_filename=          prompt(info, "OUTPUT_FERMI_PREFIX", "fermi");

  box[0]=               (int) val(prompt(info, "T", "16"));
  box[1]=               (int) val(prompt(info, "X", "4"));
  box[2]=               (int) val(prompt(info, "Y", "4"));
  box[3]=               (int) val(prompt(info, "Z", "4"));

  string action_type=             prompt(info, "ACTION_TYPE", "CLOVER");
  string inverter=                prompt(info, "INVERTER", "MINRES");
  mdp_real absolute_precision=val(prompt(info, "ABSOLUTE_PRECISION", "0"));
  mdp_real relative_precision=val(prompt(info, "RELATIVE_PRECISION", "1e-6"));
  int source_spin=      (int) val(prompt(info, "SOURCE_SPIN", "0"));
  int source_color=     (int) val(prompt(info, "SOURCE_COLOR", "0"));

  coeff["kappa_s"]=           val(prompt(info, "KAPPA_S", "0.13"));
  coeff["kappa_t"]=           val(prompt(info, "KAPPA_T", "0.13"));
  coeff["r_s"]=               val(prompt(info, "R_S", "1.0"));
  coeff["r_t"]=               val(prompt(info, "R_T", "1.0"));
  coeff["c_{sw}"]=            val(prompt(info, "C_SW", "0.0"));
  coeff["c_E"]=               val(prompt(info, "C_E", "1.0"));
  coeff["c_B"]=               val(prompt(info, "C_B", "1.0"));
  coeff["alpha_1"]=           val(prompt(info, "ALPHA_1", "0.0"));
  coeff["alpha_2"]=           val(prompt(info, "ALPHA_2", "0.0"));
  coeff["alpha_3"]=           val(prompt(info, "ALPHA_3", "0.0"));
  coeff["alpha_4"]=           val(prompt(info, "ALPHA_4", "0.0"));
  coeff["alpha_5"]=           val(prompt(info, "ALPHA_5", "0.0"));
  coeff["alpha_6"]=           val(prompt(info, "ALPHA_6", "0.0"));
  coeff["alpha_7"]=           val(prompt(info, "ALPHA_7", "0.0"));
  coeff["alpha_8"]=           val(prompt(info, "ALPHA_8", "0.0"));
  coeff["alpha_9"]=           val(prompt(info, "ALPHA_9", "0.0"));

  mpi << coeff;
		 
  int i,a;
  mdp_lattice lattice(4,box,default_partitioning0,torus_topology,0,2,false);
  gauge_field U(lattice,3);
  fermi_field psi(lattice,3);
  fermi_field chi(lattice,3);
  fermi_field phi(lattice,3);
  site x(lattice);

  if(gauge_filename=="HOT") set_hot(U);
  else if(gauge_filename=="COLD") set_cold(U);
  else U.load(gauge_filename);

  if(coeff.has_key("c_{sw}")) compute_em_field(U);

  psi=0;
  if(on_which_process(lattice,0,0,0,0)==ME) {
    cout << " This is " << ME << endl;
    x.set(0); 
    psi(x,source_spin,source_color)=1;
  }
  psi.update();
  cout << "done\n";
  
  if(action_type=="CLOVER") 
    default_fermi_action=FermiCloverActionFast::mul_Q;
#if defined(SSE2)
  else if(action_type=="CLOVERSSE") 
    default_fermi_action=FermiCloverActionSSE2::mul_Q;
  else if(action_type=="FERMILABSSE") 
    default_fermi_action=FermiFermilabActionSSE2::mul_Q;
#endif
  else
    error("unknown action");

  if(inverter=="MINRES") 
    default_fermi_inverter=MinimumResidueInverter<fermi_field,gauge_field>;
  else if(inverter=="BICGSTAB") 
    default_fermi_inverter=BiConjugateGradientStabilizedInverter<fermi_field,gauge_field>;
  else
    error("unknown inverter");

  mpi << "HERE\n";

  inversion_stats stats=mul_invQ(chi,psi,U,coeff,absolute_precision,relative_precision);

  mpi << "checking inversion...";
  chi.update();
  mul_Q(phi,chi,U,coeff);
  check_differences(psi,phi);

  mpi << "Time per mul_Q step per site=" 
      << stats.time / stats.mul_Q_steps * 1.0e6 / lattice.global_volume() 
      << " (micro secs)\n";

  char tmp[100];
  sprintf(tmp,"_%i%i.mdp", source_spin, source_color);

  fermi_filename=fermi_filename+tmp;
  chi.save(fermi_filename);

  mpi.close_wormholes();
}
